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RESEARCH  OBJECTIVES 


This  research  focuses  on  developing  a  method  for  the  preform  design  engineering  of  ma¬ 
terial  forming  processes.  In  this  research,  a  sensitivity  analysis  method  for  preform  die  shape 
design  in  material  forming  processes  is  developed  using  the  rigid  visco- plastic  finite  element 
method.  The  preform  die  shapes  are  represented  by  cubic  B-spline  curves.  The  control  points 
or  coefficients  of  B-spline  are  used  as  the  design  variables.  The  optimization  problem  is  to 
minimize  the  zone  where  the  realized  and  desired  final  forging  shapes  do  not  coincide.  The 
sensitivities  of  the  objective  function,  nodal  coordinates,  and  nodal  velocities  with  respect 
to  the  design  variables  are  developed  in  detail.  A  procedure  for  computing  the  sensitivities 
of  history-dependent  functions  is  presented.  The  remeshing  procedure  and  the  interpola¬ 
tion/transfer  of  the  history-dependent  parameters,  such  as  effective  strain,  are  stated.  The 
procedures  of  sensitivity  analysis  based  preform  die  design  are  also  described.  In  addition, 
a  method  for  the  adjustnient  of  the  volume  loss  resulting  from  the  finite  element  analysis  is 
given  in  order  to  make  the  workpiece  volume  consistent  in  each  optimization  iteration.  The 
method  developed  in  this  report  is  used  to  design  the  preform  die  shape  of  H-shaped  forging 
processes,  including  plane  strain  and  ajdsymmetric  deformations.  The  results  show  that  a 
flashless  forging  with  a  complete  die  fill  is  realized  using  the  optimized  preform  die  shape. 
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INTRODUCTION 


In  metal  forming  processes,  net-shape  manufacturing  depends  on  the  exactness  of  the  die 
design  to  a  large  extent,  especially  the  preform  die  shape  design.  Preform  die  shape  design 
engineering  includes  preform  die  and  process  sequence  designs.  They  are  the  most  impor¬ 
tant  steps  for  the  quality  control  of  products,  material  savings,  and  reduction  of  the  die 
design/manufacture  cost.  Engineering  experience  and  intuition  have  been  the  primary  tools 
in  the  forging  industry  for  determining  the  number  of  stages  and  the  preform  die  shapes. 
Meted  forming  simulations  using  finite  eleihents  have  brought  opportunity  for  development 
of  new  methods  for  process  sequence  designs. 

Kobayashi  et  al.  [1]  introduced  a  finite  element  method  that  simulated  the  metal  forming 
processes  in  reverse.  The  process  known  as  backward  tracing  was  applied  to  several  practi¬ 
cal  forging  problems  [2,3,4].  Han  et  al.  [5]  applied  an  optimization  method  in  back-tracing. 
Their  method  determined  which  nodes  to  detach  from  the  dies  based  on  nodal  velocities 
resulting  from  minimizing  the  difference  between  the  maximum  effective  strain  rate  and  the 
average  effective  strain  rate  for  the  entire  billet  domain.  The  optimization  sequence  is  im¬ 
plemented  in  each  incremental  step.  Zhao  et  al.  [6,7]  established  a  detachment  criterion  for 
backward  simulation  and  the  related  preform  design  according  to  forging  shape  complexity 
control  and  applied  this  method  in  preform  design  of  axisymmetric  deformation  problems. 
Zhao  et  al.  [8]  also  gave  an  inverse  die  contact  tracking  method  for  designing  the  preform 
shape.  The  procedure  starts  with  the  forward  simulation  of  a  candidate  preform  shape  into 
the  final  forging  shape.  A  record  of  the  boundary  condition  changes  is  produced  by  identify¬ 
ing  when  a  particular  die  segment  makes  contact  with  the  workpiece  surface.  This  recorded 
boundary  condition  change  is  then  modified  according  to  the  material  flow  characteristics. 
The  modified  boundary  condition  is  finally  used  as  the  control  criterion  for  node  detachment 
during  the  inverse  deformation  simulation. 

All  the  above  methods  use  both  forward  and  backward  simulations  and  relied  on  selecting 
the  appropriate  detachment  criteria  of  boundary  nodes  during  the  backward  simulation.  In 
each  case,  the  design  objectives  were  the  preform  or  intermediate  forging  shapes  rather  than 
directly  designing  the  dies.  So  the  preform  die  shape  then  had  to  be  designed  to  produce 
the  preform  shape  of  the  workpiece. 
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Badrinarayanan  and  Zabaras  [9,  10]  developed  a  sensitivity  analysis  method  for  large  de- 
formation  of  hyperelastic  viscoplastic  solids  that  can  be  applied  to  preform  design  problems 
in  metal  forming.  Like  the  backward  tracing,  this  method  designs  the  preform  or  inter¬ 
mediate  shape  of  the  workpiece  instead  of  preform  dies.  Their  method  was  applied  to  an 
axisymmetric  disk  upsetting  problem  where  the  preform  is  designed  such  that ,  a  final  forging 
with  a  minimum  barreling  effect  is  achieved.  The  desired  results  were  achieved,  however  the 
axisymmetric  preform  shape  had  a  concave  lateral  shape  that  is  very  difficult-to-forge. 

Fourment  et  al.  [11]  described  a  method  to  design  the  preform  tools  and  the  preform 
shapes.  The  distance  between  the  achieved  and  required  part  is  used  as  the  objective  func¬ 
tion  to  be  minimized.  The  shapes  are  discretized  using  spline  functions.  The  design  variables 
of  the  optimization  problem  are  the  displacements  of  the  selected  characteristic  points  of  the 
spline  in  the  normal  direction.  The  gradients  are  calculated  analytically  where  the  friction  on 
the  tool-workpiece  interface  is  considered  as  an  exponential  function  of  the  sliding  velocity. 
Shape  optimization  for  both  one  and  two  step  forging  operations  was  performed  using  this 
method. 

This  project  focuses  on  optimal  design  of  the  preform  die  shapes  instead  of  the  preform 
shapes.  An  optimization  method  for  preform  die  shape  design  in  metal  forming  is  developed 
using  forward  simulation  only.  The  preform  die  shapes  are  represented  using  piecewise  cubic 
B-splines  function.  The  objective  is  to  reduce  the  zone  where  the  actual  final  forging  shape 
and  desired  final  forging  shape  do  not  coincide.  B- spline  coefficients  are  considered  as  the 
design  variables  for  the  sensitivity  analysis  and  optimization  problem.  After  completing  the 
sensitivity  analysis,  the  optimization  step  is  performed  at  the  end  of  simulation  for  the  select¬ 
ed  spline  coefficients.  The  updated  spline  coefficients  define  the  new  preform  die  shapes  and 
a  complete  forging  simulation  is  again  carried  out  to  evaluate  the  sensitivity  and  objective 
function.  This  preform  die  shape  updating  procedure  is  continued  until  no  further  improve¬ 
ment  is  achieved.  By  minimizing  the  objective  function,  the  net-shape  manufacturing  of  the 
final  forging  can  be  realized  using  the  designed  preform  die  shapes.  The  required  sensitivity 
analysis  for  rigid-plastic  or  visco-rigid  plastic  deformation  problems  is  developed  in  detail, 
including  the  velocity  boundary  conditions  for  the  contact  problem.  The  friction  boundary 
condition  at  the  die- workpiece  interface  for  the  sensitivity  analysis  is  modeled  by  an  arctan¬ 
gent  function  [1].  The  remeshing  procedures  and  volume  loss  adjustment  are  also  described 
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and  examples  are  presented  using  the  sensitivity  analysis  and  optimization  method. 
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FINITE  ELEMENT  EQUATIONS 


In  the  following,  the  governing  equilibrium  equations  are  presented  for  establishing  the 
objective  function  gradients.  For  the  complete  visco-plastic  analysis,  reference  [1]  is  suggest¬ 
ed. 

The  element  strain  rate  vector  for  2D  analysis  problems  can  be  written  as: 


/  4  \ 

dvx/dx  \ 

4  _ 

dvyfdy 

4  - 

0 

\7xy  / 

\  dvx/dy  +  dvy  jdx) 

/  dvrjdr 
dvz  /  dz 

Vrlr 

\  dvz  /dr  +  dvr  jdz  > 


for  plane  strain 


for  axisymmetric  deformation 


(la) 


(16) 


Substituting  the  elemental  admissible  velocity  field  of  the  4-node  linear  element  in  equations 
(la)  and  (lb),  the  element  strain  rate  vector  is  represented  in  a  unified  form,  as 


e  = 


(ii\ 

/ 

£2 

1 

\€4/ 

1 

!  «(")  J.  TT 
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where  (fi“\t'2'*^)  velocity  components  of  the  ath  node.  They  correspond  to  Vx  and  Vy, 

respectively.  The  summation  is  over  all  four  nodes.  Fa  is  zero  for  a  plane  strain  problem  and 
the  row  of  is  is  deleted  for  plane  stress  deformation.  For  an  axisymmetric  problem,  vi  and 
represent  Vr  and  v^,  respectively,  and  Pa  becomes  qa/r,  where  r  is  the  radial  coordinate  and 
expressed  as; 


*■  =  53 

Q=1 

Strain  rates  in  the  finite  element  formulation  are  typically  computed  from 

e  =  BV 

where  B  is  the  element  strain  rate  matrix  and  V  is  the  nodal  velocity  vector  of  the  element. 
B  matrix  is  expressed  as: 


B  = 


El 

0 

Es 

0 

Es 

0 

E4 

0 

Hi 

0 

Hi 

0 

0 

Pi 

0 

P2 

0 

Ps 

0 

P4 

Pi 

El 

Hi 

Ei 

Hz 

P4 

(4) 


5 


The  elemental  stiffness  equation  of  the  forging  problem  using  FEM  can  be  expressed  as: 


K(V,X)V  +  F(V,X)-0  (5) 


where 


Kij{Y,X}  = 

/  ZPijdV  +  Q  f  CiCjdV 

(6a) 

Jv  £  Jv 

F.(V,X)  = 

[  mAj--g,ta7i“^(^— ^)<iS 

Isa  ^  ^0 

(66) 

where  K  is  the  material  and  process  dependent  nonlinear  stiffness  matrix,  F  is  the  applied 
nodal  point  force  vector.  X  is  the  nodal  coordinate  vector  of  the  element.  Q  is  a  large  positive 
constant  which  penalizes  the  dilatational  strain  in  the  rigid-plastic  formulation,  qj  are  the 
element  shape  functions  expressed  in  the  natural  coordinate  system  (£,7?).  a  is  the  effective 
stress,  a  =  a{e)  and  a  =  a‘[e/e)  for  rigid-plastic  and  rigid- viscoplastic  materials,  respectively. 
k  is  the  effective  strain  rate.  The  friction  on  the  interface  between  the  workpiece  and  dies  is 
modelled  using  the  arctangent  function  of  the  relative  sliding  velocity,  k  is  the  shear  yield 
stress  and  m  is  the  constant  friction  factor.  Ugj  is  the  relative  sliding  velocity  at  node  j  on 
the  die-workpiece  contact  interface.  The  sliding  velocity  Ug  is  approximated  using  a  linear 
shape  function: 

2 

i=i 

uq  is  a  small  positive  number  compared  to  Ug.  The  components  C,  of  matrix  P  and  vector 
C  are  constructed  by  using  the  strain  rate  matrix  B  as  follows: 


where  the  diagonal  matrix  D  is: 


Pij  —  ^ki^kn^nj 
Ci  =  Bii  -h  B2i  +  Bsi 

Cj  =  Bij  “h  B2j  +  B^j 


D  = 


(7a) 

(76) 

(7c) 
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In  metal  forming,  the  stiffness  equation  (5)  is  nonlinear.  Evaluating  stiffness  matrices  at 
the  elemental  level  from  equations  (5)  and  (6),  assembling  them  for  the  whole  workpiece,  we 
obtain  a  set  of  nonlinear  equations.  The  solution  is  then  obtained  by  using  Newton-Raphson 
method  which  consists  of  linearization  and  apphcation  of  convergence  criteria. 
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OBJECTIVE  FUNCTION  AND  DESIGN  VARIABLES 


For  a  two-dimensional  metal  forming  processes,  suppose  that  Go  and  G  represent  the 
desired  shape  and  the  actual  shape  of  the  final  forging,  respectively.  The  difference  between 
Go  and  G  is  a  function  of  the  preform  shape  and  reflects  the  exactness  of  the  preform  design. 
It  is  our  design  objective  to  make  shape  G  close  to  Go  by  redesigning  the  preform  dies.  The 
objective  function  or  difference  of  two  shapes  can  be  expressed  as  the  area  of  the  zone  where 
the  two  shapes  do  not  coincide  (see  figure  1). 

For  the  discretized  boundaries  of  the  workpiece,  suppose  there  are  N  boundary  nodes 
around  the  achieved  final  shape  G  with  coordinates  (x,,  j/j)  (i  ~  1,2,  Similarly,  the 

boundary  of  the  desired  shape  can  be  discretized  by  extending  a  line  on  the  normal  direction 
from  each  node  on  the  material  boundary  to  intersect  the  desired  shape  boundary.  This 
provides  a  second  set  of  node  coordinates  for  the  desired  shape  (xoj,yo,)  (i  =  l,2,...,jV).  By 
connecting  two  consecutive  nodes  from  each  boundary,  a  quadrilateral  element  is  formed. 
The  four  nodes  of  each  element  are  locally  numbered  in  counterclockwise  direction  as  shown 
in  figure  2(a)  and  so  the  area  of  the  element  can  be  calculated  as: 

1  1 

—  ®12j/32)  +  -(®14j/34  “  *342/14) 

where  Xij  —  Xi  Xj ,  —  y,*  yj . 

Occasionally,  the  actual  final  forging  boundary  will  intersect  the  desired  final  forging 
boundary  when  this  intersecting  area  occurs  between  two  nodes.  The  element  wiU  appear 
as  shown  in  figure  2(b).  For  this  particular  element,  its  area  is  calculated  as  follows: 

Aj  =  “(X35y25  —  *252/35)  +  “(^142/54  “  *542/14) 

where  (*5,2/5)  is  the  intersecting  point. 

We  use  the  sum  of  the  square  of  the  subarea  Aj  as  the  objective  function  V'* 

=  (8) 

i=i 

When  tjj  approaches  zero,  the  achieved  shape  G  will  be  consistent  with  the  desired  shape  Gq. 
Therefore,  the  optimization  problem  is  to  define  a  preforming  operation  that  will  minimize 
the  objective  function  V"* 
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The  shapes  of  the  preform  dies  for  a  two-dimensional  model  are  represented  using  cubic 
B-spline  functions.  The  B-spline  shapes  are  controlled  by  varying  the  coefficients  or  the 
coordinates  of  the  control  points.  For  each  control  point,  there  are  two  degrees  of  freedom 
{PziiPyi)  *  =  1,2, ....if  for  a  total  of  2K  design  variables.  For  this  unconstrained  problem,  the 
BFGS  algorithm  [12]  is  used  to  minimize  the  objective  function  ^  with  respect  to  the  design 
variables  pj. 
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SENSITIVITY  ANALYSIS 


From  the  objective  function  given  by  equation  (8),  the  gradient  of  the  objective  function 
^  with  respect  to  the  design  variable  pi  is  obtained  as  follows: 


^  =  ('=u . .i.) 

Differentiating  the  cost  function  V*  with  respect  to  the  coordinates  z,-  and  y.-  gives: 


^  dip  dxi  ^  d‘^  dyi 


(9) 


JV 


For  a  regular  element  j: 


For  an  intersecting  element  j: 
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-  =  -*35 

dy2  2 


After  convergence  of  the  finite  element  analysis,  the  velocity  field  at  the  incremental 
simulation  step  is  used  to  update  the  nodal  coordinates  using  the  following  equation: 


X(‘+At)  _  x(*)  +  v(‘^At  (10) 

where  is  the  nodal  coordinate  vector  at  time  t+ At.  is  the  nodal  coordinate  vector 

at  time  i  and  is  the  nodal  velocity  vector  at  time  t.  Differentiation  of  equation  (10)  with 
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respect  to  the  design  variables  pj  results  in 


+  — At 


dpi 


dpi  dpi 


(11) 


At  the  initial  state  (t  =  o): 


dX(°) 

dpi 


=  0 


(12) 


It  can  be  seen  from  equations  (9)  and  (11)  that  the  gradients  of  the  objective  function 
with  respect  to  the  design  variables  can  be  computed  once  the  sensitivities  of  the  nodal 
velocities  with  respect  to  the  design  variables  are  available. 


DETERMINATION  OF 


dpi 


The  nodal  velocity  sensitivities  V,p,  at  time  f  are  obtained  by  differentiating  the  equilib¬ 
rium  equation  (5)  with  respect  to  the  design  variable  pi  at  the  elemental  level. 


,,  dT^dX  ,dKdX^^  dFdX, 

'^ax^dpi~  ^  ax  dpi^  ^  dX  dpi^ 


(13a) 


where  X  and  V  are  the  nodal  coordinate  and  velocity  vectors  respectively  for  the  element  at 
time  i. 

Equation  (13a)  can  be  simply  rewritten  as 


—  ®'>pi 


(136) 


where  R  is  the  elemental  stiffness  matrix  sensitivity  (8  x  8),  F,p,  is  the  elemental  force  vec¬ 
tor  sensitivity  and  V.p,  is  the  sensitivity  vector  of  the  nodal  velocities  with  respect  to  the 
Ith  design  variable.  The  components  of  the  matrix  R  and  vector  F,p,  are  expressed  as  follows: 
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=  E  E  /  (ill  -  ^)\PinVnVmPmidV  +  Kij  +  ^  {i,  j  =  1.  2,  ...,8)  (14a) 

Ti  =  l  mnl  ^  ^ 


=  -(E  +  E  E  S^Sr^i)  (*  =  1. 2.  8) 


(146) 


^idx„api  f^.^^dxnapi 

After  evaluating  the  sensitivities  of  the  stiffness  matrix  and  the  nodal  force  vector  at  the 
elemental  level  using  equations  (13)  and  (14),  they  are  assembled  for  the  whole  workpiece. 
A  set  of  simultaneous  linear  algebraic  equations  is  obtained  as 


For  a  given  initial  guess  of  the  preform  die  shape,  the  finite  element  analysis  of  the 
preform  stage  is  performed.  In  each  incremental  simulation  step,  the  sensitivities  of  the 
node  velocities  to  the  design  variables  are  obtained  by  solving  equation  (15)  after  the  finite 
element  solution  in  this  step  has  converged.  The  sensitivities  of  the  nodal  coordinates  with 
respect  to  the  design  variables  are  dependent  on  the  deformation  history  and  so  are  updated 
using  equation  (11).  Once  the  simulation  of  the  preform  stage  is  finished,  the  simulation  of 
the  final  forging  stage  is  started.  During  the  simulation  of  the  final  stage,  the  sensitivities  of 
the  boundary  node  velocities  with  respect  to  the  design  variables  are  zero  since  the  final  die 
shape  does  not  change.  When  the  final  stage  simulation  is  finished,  the  resulting  sensitivi¬ 
ties  of  the  nodal  coordinates  with  respect  to  the  design  variables  are  used  to  calculate  the 
objective  function  gradients  using  equation  (9)  and  the  objective  function  is  calculated  using 
equation  (8).  After  a  complete  simulation  including  the  preform  and  final  stages  is  finished, 
the  optimization  program  DOT  [13]  is  called  to  provide  new  control  point  coordinates  for 
the  preform  dies.  Usually,  no  upper  and  lower  bounds  are  imposed  on  the  design  variables 
except  in  particular  cases  such  as  when  a  control  point  of  the  preform  die  must  be  fixed  if 
half  or  quarter  symmetry  is  used  in  the  analysis.  At  this  point,  the  optimization  program 
also  checks  to  see  if  the  optimality  conditions  are  satisfied.  If  the  optimality  conditions  are 
not  satisfied,  the  preform  die  shapes  are  updated  using  the  new  control  points  and  the  next 
optimization  step  begins  with  the  initial  state.  If  the  optimality  conditions  are  satisfied  with 
the  current  shapes,  the  design  objectives  are  met.  The  design  optimization  flow  chart  for 
the  two-stage  forging  process  is  shown  in  figure  3. 


DETERMINATION  OF 


dKii  dFi  dFi 
dXn  ^  dXn^  dVj 


A,ccording  to  equations  (6a)  and  (6b),  the  derivatives  and  in  equation  (14) 

at  the  elemental  level  are  developed  as  follows: 

=  /  mk-qiqj  2  I (16) 
dvj  Js  IT  %  +  {qkU,ky 

I?  =  EE  /  -^^M^Vktan-H^)dS+  f  mk-qMn-\^)^  (17) 

°  ^  ifei •'5  de  e  dxn  tto  isc  ^  '  ito 


f.j  f  Add-  ff-,1  a(v^pv)„  r  ^  f  dCi^  ^  f  ^  sc, 

—  =  /  (t^-^  ^-="-5 - -PijdV+  /  ^^dV  +  Q  /  ^CjdV  +  Q  /  Ci-^c 

Kn  Jy  £  oe  ^  2e  OXn  Jy  £  OXn  Jy  dXn  Jy 

M  _  £\.  I  Vv  /  «vv\ 


Jy  £  dXn  Jy 
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^y^here  1^,  1^,  ^  and  can  be  obtained  according  to  equations  (7a)  through  (7c). 

aXn  ^  OXn  '  O^n 

^  _  dBu  ^  dB2i  ^  dB^i 

dxfj  d^Yi 


dCj  ^  dBij  ^  dB2j  ^  dBsj 

dXn  dXn  dXn  ^®n 


(196) 


dPij 

dXn 


1 _ 1  -*.  —  1  ^ 


dB, 


rnj 


dXn 


g(V'^PV) 

dx„ 


8  8 


=  EE». 


^P mk 
dXn 


Vk 


(19c) 

(19d) 


The  derivative  matrix  can  be  obtained  by  calculating  from  the  B 

matrix  given  by  equation  (4).  The  components  of  are  expressed  as  follows: 
dEildxj  =  -EiEj  (i,j  =  1,2,3,4) 

dHi/dyj^-HiHj  (i,i=  1,2,3,4) 


dEi/dyj  =Aij/8\3\-  EiHj  (t.j  =  1,2,3,4) 
dHi/dxj  =  ->l.j/8|J|  -  HiEj  {i,  j  =  1, 2, 3, 4) 
dPi/dxj  =  {hj—  1)2)3, 4) 

dPi/dyj=0  (i,i=l,2,3,4) 

where  |J|  is  the  determinant  of  Jacobian  matrix  and  expressed  as: 


l^i  =  X  [(*133/24  -  *243/13)  +  (*343/12  -  *123/34 )^  +  (*233/14  “  *143/23)»7]  (20) 

o 

where  x.^  =  x,-  -Xj, y^  =  yi-yj  {i,j  =  1, 2, 3,4).  For  a  plane  strain  problem,  dPi/dxj  =  dPi/dyj  =  0 
{i,j  =  1,2, 3, 4).  The  coefficients  Aij  are  given  by  a  matrix  A 


/  0  ij-1  (-ri  l-(  \ 

I  1-jj  0  -(1+0  i  +  V  I 

[fj-i  1+e  0  -(1+1?) 

V^-I  -(^  +  1?)  1  +  1?  0  / 


Interestingly,  the  coefficient  matrix  A  is  an  antisymmetric  matrix  because  A’’  =  -A. 

It  should  be  noted  that  the  differential  volume  dV  and  area  dS  are  also  dependant  on  the 
nodal  coordinate  vector  X.  By  using  the  natural  coordinates  (-1  <  ^  <  1,  -1  <  i?  <  l)  in  the 
two-dimensional  space,  and  can  be  represented  as  follows. 

For  a  plane  strain  problem: 


d(dS)  ^  g|J,| 

dXji  dXn 
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a(dv)  _  s|j| 

_ —  =  _ — d^df] 

uXfi 

where  |Js|  is  the  determinant  of  the  Jacobian  of  the  coordinate  transformation  matrix  on  the 
die- workpiece  interface. 


|Js|  =  +  (yi  -  y2y  (21) 

where  (xi,yi)  and  (x2,y2).  are  the  coordinates  of  the  two  nodes  of  the  element  in  contact  with 
the  die  surface. 

For  an  axisymmetric  problem: 


djdS) 

dx„ 


where  r  is  the  radial  position  of  the  integration  point  and  defined  by  equation  (3) 


dr 

dX2i-l 


dr 

dX2i 


=  0 


(i=l,2,3,4) 


and  can  be  obtained  by  using  equations  (20)  and  (21),  respectively. 


BOUNDARY  CONDITIONS 

On  the  friction  boundary,  the  traction  is  prescribed  in  the  tangential  direction  and  the 
velocity  is  prescribed  in  the  normal  direction  to  the  surface.  When  the  interface  direction 
is  inclined  with  respect  to  the  global  coordinate  axis,  the  coordinate  transformation  of  the 
stiffness  matrix  upon  the  inclined  direction  is  necessary  in  order  to  impose  the  mixed  bound¬ 
ary  conditions.  The  local  coordinate  system  is  defined  as  the  inclined  boundary  coordinate 
system  (see  figure  4).  The  velocity  boundary  conditions  of  the  ith  node  in  contact  with  the 
dies  are: 

=  VL.n=  (22) 

where  is  the  velocity  component  of  the  node  *  in  the  normal  direction  of  the  interface 
surface.  V*e  is  the  die  velocity  vector,  n  is  the  unit  normal  on  the  interface  surface.  /3 
is  measured  from  the  x  axis  in  the  global  coordinate  system  to  the  x'  axis  of  the  local 
coordinate  system  in  counterclockwise  direction.  For  a  B-spline  function  defined  by  y  =  y(p/,  x) 
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(/  =  1, 2, 2K),  the  slope  of  the  B-spline  curve  is  ^  Thus,  the  following  relations 

are  obtained 


sinfS  = 


Vx 

Vi  +  Vx 


C05/3  = 


+  3/1 

Differentiating  equation  (22)  with  respect  to  the  design  variable  pi  gives  the  following 
relationship: 

(23) 


^  yT 


dpi 


dsinfi  dcosfi 


^liere  ^  is  the  sensitivity  of  the  normal  of  the  die  surface  to  the  design  variable  p/.  It  can 
be  obtained  by  differentiating  the  equation  which  defines  the  normal  of  the  B-spline  curve. 


where  y^p,  =  If  the  node  i  is  not  in  contact  with  part  of  the  dies  that  is  being  optimized, 
or  if  the  dies  are  not  optimized  (for  instance,  during  the  final  stage  forging),  then  ^  =  0. 

Usually,  is  equal  to  zero.  So  equation  (23)  can  be  simplified  into 

dcosP  i  ^ 

-3 —  =Vdiey—E - 

dpi  *  dpi 

The  above  boundary  conditions  should  be  imposed  in  the  local  coordinate  system  and 
the  sensitivity  equations  (15)  are  solved  in  the  local  coordinate  system.  When  the  solution 
is  obtained,  it  is  transformed  back  to  the  global  coordinate  system. 

CUBIC  B-SPLINE  AND  ITS  DIFFERENTIATION 

Given  the  control  points  (?*<,?»;)  (»  =  1, 2, ....  JC),  the  design  variables  are  p  =  (pi,  P2, P2jc)  = 
(Pxi,Py,,",PxK>P»K)-  The  (®,y)  coordinates  of  the  ith  piecewise  cubic  B-spline  function  can  be 
expressed  by  a  parametric  function  as  [14] 


a  =  1  [(1  -  s)3p^ .  +  (4  -  Gs^  +  3«3)p^.^^  +  (1  +  3s  +  3s''  -  3s^)px,^,  +  s^Px.+a]  (26o) 

6 

y  =  ■r[(l  ~  *)^P»,-  -l-  (4  —  6s^  +  3s^)pyj^,  +  (1  -|-  3s  +  3s^  —  3s®)pyi^2  -t-  s  Pyi^.3]  (266) 

6 

where  0  <  s  <  1. 

Differentiation  of  the  B-sphne,  y*  can  be  obtained  using  the  following  equation: 


_  ^  _  dy  .dx 
dx  ds  ds 


(27) 
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where 


I;  =  ^[-(1  -  ofvvi  +  (-4s  +  3s2)pa;^,  +  (1  +  2s  -  +  s2pj,,+3]  (28o) 

=  ^[-(1  -  +  (-4s  +  ^s^)Pxi+,  +  (1  +  2s  -  .3s2)p3:,+,  +  S Vi+al  (286) 

The  differentiations  of  the  slope  of  B-spline  to  the  design  variables  or  control  points  (pr.,Pyi) 
are  as  follows; 


d^y  _ 

(l-«)2 

dx  dpy- 

9  dx 
^  ds 

d^y 

—45  -h  35^ 

9xdpy.^, 

9  dx 

d^y 

1  -1-  2s  -  3s2 

dxdpy-,^,^ 

t\dx 

^ds 

d^y 

_  s2 

9xdpy,^,  2|f 

d'^y  __ 

(1  -  s)2  dy 

2{%f  ds 

d^y  _ 

—45  +  35^  dy 

2{^y  ds 

d^y  _ 

1  -h  2s  -  3s2  dy 

dxdpxi^x 

2(f)'  5s 

(29a) 

(296) 

(29c) 

(29rf) 

(30a) 

(306) 

(30c) 


d^y  _  dy 

5a9pxi+3  2(|f  )2  as 


(sod) 


Substituting  yx  (equations  (27)  and  (28))  and  (equations  (29)  and  (30))  in  equation 
(24),  is  solved.  Then  substituting  equation  (24)  in  equation  (23)  or  (25)  gives  the  velocity 
boundary  conditions  for  the  sensitivity  analysis. 


REMESHING 

In  practical  forming  processes,  large  deformations  eventually  lead  to  indeterminant  ele¬ 
ments  when  the  determinant  of  Jacobian  matrix  becomes  negative.  Therefore,  a  new  mesh 
of  the  workpiece  must  be  defined  and  the  history-dependant  variables  must  be  transferred 
to  the  new  mesh  system.  The  history-dependant  field  variables  are  effective  strains,  node 
temperatures  (only  for  non-isothermal  analysis)  and  the  sensitivities  of  the  node  coordinates 
with  respect  to  the  design  variables.  These  values  must  be  defined  on  the  new  mesh  by 
interpolation. 

Sensitivities  of  the  node  coordinates  with  respect  to  the  design  variables  are  given  at  the 
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node  points.  Thus,  it  is  assumed  that  the  distributions  of  node  coordinate  sensitivities  within 
the  workpiece  domain  can  be  expressed  by  using  the  element  shape  functions.  Interpolation 
is  done  by  evaluating  the  sensitivities  at  the  new  node  locations. 


ADJUSTMENT  OF  VOLUME  LOSS 

A  small  volume  loss  of  the  workpiece  due  to  geometry  update  within  a  finite  time- 
increment  is  inevitable.  In  addition,  the  amount  of  volume  loss  will  vary  due  to  remeshing. 
Limiting  the  volume  loss  within  a  small  percentage  of  the  total  deforming  volume  is  a  major 
consideration  in  the  prediction  of  proper  die  fiU  and  defect  formation,  which  are  important 
in  process  design.  In  a  conventional  finite  element  simulation,  the  amount  of  volume  loss  is 
controUed  by  hmiting  the  maximum  allowable  time-increment. 

The  initial  billet  volume  is  determined  to  be  equal  to  the  final  forging  volume.  However, 
with  volume  loss,  the  achieved  final  forging  volume  V  is  always  smaller  than  the  initial  billet 
volume  Vo-  The  objective  function  depends  on  the  area  of  the  zone  where  two  final  forging 
shapes  do  not  coincide.  Therefore,  the  objective  function  is  affected  adversely  by  the  volume 
loss.  More  importantly,  the  magnitude  of  volume  loss  may  vary  from  one  optimization 
iteration  to  another.  This  results  in  convergence  problems  of  the  optimization  algorithm. 

To  support  design  optimization,  a  volume  loss  adjustment  procedure  is  incorporated  in 
the  simulation  process.  After  each  time-increment,  the  y-coordinates  of  the  nodes  in  contact 
with  the  die  surfaces  are  adjusted  to  ensure  the  workpiece  volume  equal  to  the  initial  volume 
Vq.  At  the  same  time,  the  die  position  is  also  adjusted  by  the  same  distance. 

Figure  5  shows  an  element  in  contact  with  the  upper  die.  Suppose  that  the  coordinates 
of  node  4  are  (*4,2/4)  and  coordinate  adjustment  amount  in  y-direction  is  Ay.  The  new 
coordinates  of  node  4  are  (*4,  y4+Ay).  For  axisymmetric  problems,  the  volume  of  the  element 
1234  is: 

Fe  =  -Ci(*322/i2  —  ®12lfe2)  +  -C'2(*14l/34  —  »34S>14) 

where  Ci  =  ^(*1  -f  *2  -I-  *3),  =  ^(*1  -I-  *3  -I-  *4).  For  plane  strain  problems,  Ci  =  C'2  =  1. 

The  volume  of  the  element  1234'  is: 

Ve  =  2^l(®32l/l2  —  »12!/32)  +  2^2(»14l/34'  —  »34!/l4') 

The  volume  increase  due  to  the  adjustment  of  the  node  4  is: 

AV4  =  V;  -  F  =  ic'2(*34  -  *14)  •  Ay  =  54  •  Ay 
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By  accumulating  the  volume  increase  due  to  the  nodal  coordinate  adjustment  over  all  the 
nodes  in  contact  with  the  dies,  the  total  volume  increase  AFt  is  obtained. 

AVT  =  Ay^Si 

Letting  the  total  volume  increase  due  to  the  nodal  coordinate  adjustment  equal  to  the 
volume  loss  (Fq  -  resulting  from  finite  element  analysis,  we  obtain  a  factor  for  adjusting 
the  y-coordinates  of  the  nodes  in  contact  with  the  dies. 

,  AFt  Fo  -  F 

^~ESi~  ESi 

Numerical  analysis  examples  show  that  the  adjustment  amount  Ay  for  volume  loss  in 
each  time- increment  step  is  very  small.  Usually  the  ratio  of  Ay  to  the  time-increment  is 
about  0.001  ~  0.0001.  Such  a  small  adjustment  does  not  affect  the  material  flow  patterns.  But 
after  the  adjustment,  the  volume  constancy  is  realized  in  every  time-increment  step.  The 
final  forging  volume  is  exactly  equal  to  the  initial  volume  in  various  optimization  iterations 
and  the  convergence  of  optimization  is  improved  significantly. 
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DESIGN  EXAMPLES 


The  sensitivity  analysis  presented  here  are  used  to  design  the  preform  die  shape  of  H- 
shaped  forging  in  both  plane  strain  and  axisymmetric  deformation  modes.  The  goal  is  to 
design  the  preform  die  shape  such  that,  after  the  final  forging  stage,  a  flashless  forging  with 
a  complete  die  fiU  is  obtained. 

In  each  case,  the  forging  process  is  isothermal.  A  non-strain  hardening  material  having 
the  constitutive  relation  a  =  with  lo  =  33.99Mpo  was  used  in  this  example.  A  con¬ 

stant  shear  friction  factor  was  assumed  between  the  workpiece  and  the  dies.  The  lower  die 
is  stationary  and  the  upper  die  velocity  is  -l.Omm/s.  Since  the  forging  is  symmetric  about 
the  vertical  and  horizontal  axes,  only  one  quarter  of  the  cross  section  is  considered  for  the 
finite  element  analysis.  The  selection  of  the  total  number  of  control  points  representing  the 
preform  die  shapes  depends  on  the  complexity  of  the  final  die  shapes.  If  the  final  die  shapes 
are  complex,  more  control  points  should  be  used.  However,  additional  control  points  need 
more  computation  time.  In  this  example,  ten  control  points  {pxuPyi)  (p®!  >  p®,  >  ...  >  p®io) 
were  used  to  represent  the  preform  die  shape.  Therefore,  there  are  seven  spline  pieces  with 
twenty  design  variables  p/  (i  =  1,2,  ...,20).  Phantom  control  points  are  assigned  at  the  end 
points  of  the  preform  die  p^^  =  -p*,  =  6  and  p®,  =  0  since  the  splines  do  not  interpolate  to 
the  end  points.  Therefore,  side  constraints  are  imposed  on  design  variables  p®i,Pxj  and  px, 
during  the  optimization. 


PLANE  STRAIN  DEFORMATION 

In  this  case,  the  initial  guess  of  the  preform  die  shape  is  a  flat  die.  The  constant  shear 
friction  factor  m  is  taken  as  0.2.  The  height  of  the  preform  shape  at  vertical  symmetry  axis 
(®  =  O)  is  kept  same  for  each  optimization  iteration.  This  height  determines  the  stroke  of 
the  preforming  stage.  Figure  6  shows  the  resulting  preform  shapes  and  corresponding  final 
forging  shapes  at  various  optimization  iterations.  Figure  7  shows  the  comparison  of  the 
three  final  forging  shapes.  It  can  be  seen  that  the  optimized  shape  is  almost  same  as  the 
desired  shape.  Figure  8  shows  the  objective  function  history  over  the  optimization  process. 
Using  the  flat  die  initially  gives  a  final  forging  with  flash  and  incomplete  die  fill.  After  four 
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optimization  iterations,  the  die  cavity  is  completely  filled  and  a  flashless  final  forging  is  ob¬ 
tained.  The  objective  function  value  is  reduced  to  O.lTSmm'*  compared  to  a  starting  value  of 
17.096mm''.  For  the  fifth  iteration,  the  preform  die  has  a  relatively  deep  cavity  which  can  also 
give  a  flashless  final  forging  shape  but  the  objective  function  has  increased  slightly.  This  is 
due  to  elements  overlapping  the  curved  die  surface.  This  situation  can  be  resolved  by  using  a 
finer  mesh.  From  the  fifth  iteration  on,  the  objective  function  changes  very  little.  Therefore, 
the  preform  die  shape  at  either  the:  fourth  or  seventh  iteration  can  be  selected  as  the  optimal 
die  shape.  Figure  9  shows  how  the  preform  die  shapes  evolve  with  each  optimization  step. 


AXISYMMETRIC  DEFORMATION 

In  this  case,  the  desired  final  forging  shape  cross  section,  friction  factor  m  and  initial 
guess  of  the  preform  die  shape  are  same  as  those  in  plane  strain  deformation.  The  mesh  de¬ 
formation  patterns  at  the  end  of  the  preforming  and  final  forging  stages  at  each  optimization 
iteration  are  shown  in  figure  10.  Figure  11  shows  the  comparison  of  the  three  final  forging 
shapes.  The  optimized  shape  is  close  enough  to  the  desired  shape.  The  optimization  itera¬ 
tion  history  is  shown  in  figure  12.  Similar  to  the  previous  case,  for  the  initial  flat  preform 
die,  the  achieved  final  forging  has  a  flash  although  it  is  small  and  the  final  die  cavity  is  not 
completely  fiUed.  The  forging  shape  does  not  satisfy  the  dimensional  requirement.  After 
five  design  iterations,  the  final  die  cavity  is  completely  fiUed  and  a  flashless  final  forging 
is  obtained.  The  objective  function  is  reduced  horn  9.309mm^  to  0.210mm'*.  In  this  case, 
remeshing  was  needed.  The  preform  die  shape  iteration  history  for  axisymmetric  problem  is 
shown  in  figure  13.  From  these  two  analysis  cases,  it  can  be  concluded  that  the  numerical 
optimization  method  developed  in  this  work  is  very  effective  in  realizing  a  net-shape  forging 
process. 

The  optimization  method  and  related  sensitivity  analysis  presented  in  this  work  have 
been  incorporated  into  a  finite  element  code  developed  specifically  for  analyzing  metal  form¬ 
ing  processes.  It  should  be  emphasized  that,  in  fact,  the  experimental  verification  of  the 
preform  die  shape  design  method  developed  in  this  work  is  equivalent  to  the  verification  of 
the  finite  element  code.  Once  a  preform  die  shape  is  obtained  from  an  optimization  iteration, 
the  finite  element  analysis  is  performed  to  see  the  material  flow  and  the  workpiece  shape  at 
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the  end  of  preforming  and  final  forging  stages.  Therefore,  the  finite  element  analysis  veri¬ 
fies  the  preform  die  design  in  every  optimization  iteration  until  the  achieved  shape  is  close 
enough  to  the  desired  shape.  The  results  from  the  finite  element  code  used  in  this  report 
were  found  to  correlate  well  with  the  commercial  code  ANTARES  [15].  It  was  also  validated 
for  forging,  rolling  and  extrusion  simulations  and  experiments.  The  comparison,  validation 
and  experiment  show  that  the  finite  element  code  is  reliable.  Consequently,  the  results  of 
preform  die  shape  design  in  this  work  are  reasonable. 
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CONCLUSIONS 


An  optimization  and  sensitivity  analysis  methods  for  preform  die  shape  design  in  realiz¬ 
ing  net-shape  final  forging  processes  was  developed.  The  method  includes  the  determination 
of  the  objective  function,  sensitivity  analysis  and  the  velocity  boundary  conditions.  The 
optimization  procedures,  remeshing  and  the  adjustment  of  volume  loss  due  to  finite  element 
analysis  were  also  introduced.  The  sensitivity  of  the  objective  function  is  calculated  by  the 
accumulated  sensitivity  of  the  nodal  coordinates  to  the  design  variables  throughout  an  entire 
simulation  including  the  preforming  and  final  forging  processes.  The  method  was  applied 
to  design  the  preform  die  shape  in  H-shaped  forging  processes,  including  plane  strain  and 
ajcisymmetric  deformations.  The  preform  shapes  were  obtained  from  the  preforming  stage 
using  optimized  preform  die  shapes.  After  the  final  stage  forging,  a  flashless  forging  with 
complete  die  fill  was  produced.  The  design  results  are  consistent  with  the  desired  final 
shapes.  This  procedure  indicates  that  the  computed  sensitivity  fields  and  the  results  of  the 
preform  die  design  are  satisfactory.  Each  optimization  step  required  just  one  complete  forg¬ 
ing  simulation  and  the  procedure  is  very  efficient.  The  method  is  very  effective  in  realizing 
a  net-shape  forging.  Future  research  needs  to  include  the  optimization  of  other  important 
aspects,  such  as  uniform  deformation  and  energy  requirement,  so  as  to  realize  multi-objective 
optimization  design  of  metal  forming  processes. 

This  research  strictly  addresses  the  preform  die  shape  design.  The  initial  billet  stock 
size  (aspect  ratio)  is  assumed  empirically.  The  different  stock  size  would  result  in  different 
preform  die  shapes.  So  optimization  of  starting  stock  size  has  become  a  topic  for  future 
research.  This  could  be  incorporated  by  one  additional  design  variable,  say  height  for  ex¬ 
ample,  since  the  required  volume  is  known.  The  method  of  computing  shape  sensitivity  for 
this  design  variable  need  to  be  derived  and  implemented  in  a  manner  that  is  valid  for  both 
single  and  multi-stage  operations. 

This  research  concentrates  on  the  preform  die  shape  of  the  two-dimensional  problem- 
s.  The  method  and  ideas  can  also  be  applied  to  the  preform  die  shape  design  for  three- 
dimensional  problems,  The  objective  function  can  be  expressed  as  the  volume  of  the  zone 
where  the  desired  and  actually  achieved  shapes  do  not  coincide.  The  preform  die  shapes 
are  modeled  using  surface  splines.  The  control  points  of  the  surface  splines  are  used  as  the 
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design  variables.  For  three-dimensional  problems,  the  computer  time  may  become  the  major 
concern  due  to  a  large  number  of  elements  and  design  variables. 
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Fig  1.  .Desired  and  actually 
achieved  final  forging  shapes. 


(a)  (b) 


Fig  2.  Elements  for  calculating  the  objective  function, 
(a)  Regular  element,  (b)  Intersecting  element. 
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Fig  3.  Sensitivity  analysis  based  preform  die  shape  design 
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Fig  4.  Velocity  boundary  condition 
of  a  node  in  contact  with  die. 
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Fig  6.  Prefoim  and  final  forging  shapes  at  various  optimization 
iterations  in  plane  strain  deformation. 
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Fig  7.  Comparison  of  the  final  forging  shapes 
in  plane  strain  deformation. 


OBJECTIVE  (MM"M) 

0  B.O  8.0  12.0  16 


O 


8.  Objective  function  vs.  optimization 
iteration  in  plane  strain  deformation. 
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Fig  9.  Iteration  history  of  the  preform  die 
shapes  in  plane  strain  deformation. 
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Fig  10.  Preform  and  final  forging  shapes  at  various  optimization 
iterations  in  aodsymmetric  deformation. 
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Fig  12.  Objective  function  vs.  optimization 
iteration  in  axisymmetric  deformation. 
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Fig  13.  Iteration  history  of  the  preform  die 
shapes  in  aodsymmetric  deformation. 


